function [b,se,n,r2]=ols(x,y)

[n,k] = size(x);

X  = [ones(n,1) x];
XX = X'*X;
b = XX\(X'*y);
yhat = X*b;
ehat = y-yhat;
sig2= ehat'*ehat/(n-k-1);
vb = sig2*inv(XX);
se = sqrt(diag(vb));

ybar = mean(y);
SST = (y-ybar)'*(y-ybar);
SSR = ehat'*ehat;
r2 = 1-SSR/SST;